El avance en las tecnologías de secuenciación en los últimos años, ha permitido estudiar mayor número de secuencias genéticas provenientes de distintos organismos. Asímismo, los avances en el almacenamiento de la información secuenciada ha favorecido el análisis masivo de datos.
El sufijo “omica” refiere el estudio completo de dicho nivel molecular.
La secuenciación de transcriptomas completos permite conocer los genes que expresa una célula o un conjunto de células en un momento determinado y bajo condiciones particulares.
Algunas de las preguntas que pueden responderse con el análisis de transcriptomas son:
La técnica que permite realizar la secuenciación de los transcriptomas es la secuenciación del RNA (RNA-seq).
Una vez secuenciadas las bibliotecas se genera una cantidad masiva de datos ¿qué a hacer con toda la información obtenida?
La cantidad de datos que se generan a partir de la secuenciación de RNA es inmensa. El análisis de expresión diferencial requiere seguir el siguiente pipeline o flujo de trabajo:
pipeline
Al finalizar el proceso de estimación de la abundancia, RSEM nos devuelve los siguientes archivos por muestra:
## t5.genes.results
## t5.isoforms.results
Y tienen el siguiente contenido de archivos de texto plano:
## gene_id transcript_id(s) length effective_length expected_count TPM FPKM
## ENSG00000000003 ENST00000373020,ENST00000494424,ENST00000496771,ENST00000612152,ENST00000614008 2063.22 1892.73 369.00 19.99 14.14
## ENSG00000000005 ENST00000373031,ENST00000485971 940.50 770.06 0.00 0.00 0.00
## ENSG00000000419 ENST00000371582,ENST00000371584,ENST00000371588,ENST00000413082,ENST00000466152,ENST00000494752 1086.32 915.83 853.00 95.51 67.55
Se requieren importar los datos generados (genes.results o .transcripts.results) por RSEM a R. Para ello se requiere de la librería de tximport [(??? ctbTximport2017). Esta librería importa archivos de cuentas generados por:
Se requiere de una tabla de metadatos:
## Sample Cell Treatment Replicate unique_id
## 1 t1 A Control 1 A_Control_1
## 2 t2 A Control 2 A_Control_2
## 3 t3 A Control 3 A_Control_3
## 4 t4 M Control 1 M_Control_1
## 5 t5 M Control 2 M_Control_2
## 6 t6 M Control 3 M_Control_3
## 7 t7 A Verafinib 1 A_Verafinib_1
## 8 t8 A Verafinib 2 A_Verafinib_2
## 9 t9 A Verafinib 3 A_Verafinib_3
## 10 t10 M Verafinib 1 M_Verafinib_1
## 11 t11 M Verafinib 2 M_Verafinib_2
## 12 t12 M Verafinib 3 M_Verafinib_3
Y crear una ruta hacia la ubicación de los archivos:
##Se requiere crear una ruta al directorio en donde se encuentran los archivos
files <- file.path("../RSEM/RSEM/", samples$Sample, paste0(samples$Sample, ".genes.results"))
names(files) <- samples$unique_id
##Se corrobora que la ruta fue creada y contiene los nombres de los archivos correctamente
data.frame(files = files, Exists = file.exists(files))## files Exists
## A_Control_1 ../RSEM/RSEM//t1/t1.genes.results TRUE
## A_Control_2 ../RSEM/RSEM//t2/t2.genes.results TRUE
## A_Control_3 ../RSEM/RSEM//t3/t3.genes.results TRUE
## M_Control_1 ../RSEM/RSEM//t4/t4.genes.results TRUE
## M_Control_2 ../RSEM/RSEM//t5/t5.genes.results TRUE
## M_Control_3 ../RSEM/RSEM//t6/t6.genes.results TRUE
## A_Verafinib_1 ../RSEM/RSEM//t7/t7.genes.results TRUE
## A_Verafinib_2 ../RSEM/RSEM//t8/t8.genes.results TRUE
## A_Verafinib_3 ../RSEM/RSEM//t9/t9.genes.results TRUE
## M_Verafinib_1 ../RSEM/RSEM//t10/t10.genes.results TRUE
## M_Verafinib_2 ../RSEM/RSEM//t11/t11.genes.results TRUE
## M_Verafinib_3 ../RSEM/RSEM//t12/t12.genes.results TRUE
Cuidado: Revisar bien el orden (jerarquía) y nombre de los directorios y archivos.
Se importan los archivos empleando tximport:
## reading in files with read_tsv
## 1 2 3 4 5 6 7 8 9 10 11 12
## List of 4
## $ abundance : num [1:60675, 1:12] 16.8 0 118.2 9.8 19.4 ...
## ..- attr(*, "dimnames")=List of 2
## .. ..$ : chr [1:60675] "ENSG00000000003" "ENSG00000000005" "ENSG00000000419" "ENSG00000000457" ...
## .. ..$ : chr [1:12] "A_Control_1" "A_Control_2" "A_Control_3" "M_Control_1" ...
## $ counts : num [1:60675, 1:12] 338 0 1101 319 455 ...
## ..- attr(*, "dimnames")=List of 2
## .. ..$ : chr [1:60675] "ENSG00000000003" "ENSG00000000005" "ENSG00000000419" "ENSG00000000457" ...
## .. ..$ : chr [1:12] "A_Control_1" "A_Control_2" "A_Control_3" "M_Control_1" ...
## $ length : num [1:60675, 1:12] 1981 770 915 3201 2306 ...
## ..- attr(*, "dimnames")=List of 2
## .. ..$ : chr [1:60675] "ENSG00000000003" "ENSG00000000005" "ENSG00000000419" "ENSG00000000457" ...
## .. ..$ : chr [1:12] "A_Control_1" "A_Control_2" "A_Control_3" "M_Control_1" ...
## $ countsFromAbundance: chr "no"
## [1] "matrix"
## A_Control_1 A_Control_2 A_Control_3 M_Control_1 M_Control_2
## ENSG00000000003 338.00 297.00 336.00 492.00 369.00
## ENSG00000000005 0.00 0.00 0.00 0.00 0.00
## ENSG00000000419 1101.00 1199.00 1054.00 1122.00 853.00
## ENSG00000000457 319.26 263.52 342.49 421.71 291.48
## ENSG00000000460 454.74 609.48 416.51 502.29 380.52
## ENSG00000000938 0.00 0.00 0.00 0.00 0.00
## M_Control_3 A_Verafinib_1 A_Verafinib_2 A_Verafinib_3
## ENSG00000000003 367.00 327.00 295.00 309.00
## ENSG00000000005 0.00 0.00 0.00 0.00
## ENSG00000000419 917.00 860.00 922.00 889.00
## ENSG00000000457 281.16 361.74 372.45 379.61
## ENSG00000000460 631.84 183.26 190.55 146.39
## ENSG00000000938 1.00 1.00 2.00 1.00
## M_Verafinib_1 M_Verafinib_2 M_Verafinib_3
## ENSG00000000003 594.00 609.00 633.0
## ENSG00000000005 1.00 3.00 1.0
## ENSG00000000419 787.00 838.00 762.0
## ENSG00000000457 483.88 492.25 436.6
## ENSG00000000460 216.12 245.75 228.4
## ENSG00000000938 0.00 0.00 0.0
Para el análisis de expresión diferencial requerimos que los valores de las cuentas sean números enteros. Recodificar los valores a integer y convertir la matriz de cuentas a un data frame:
storage.mode(counts) <- "integer"
##Recodificar la matriz de cuentas a data frame
counts <- as.data.frame(counts)
class(counts)## [1] "data.frame"
## [1] "A_Control_1" "A_Control_2" "A_Control_3" "M_Control_1"
## [5] "M_Control_2" "M_Control_3" "A_Verafinib_1" "A_Verafinib_2"
## [9] "A_Verafinib_3" "M_Verafinib_1" "M_Verafinib_2" "M_Verafinib_3"
Filtrado de cuentas Filtrar los datos es muy útil. Los genes cuyos valores de cuentas son igual a 0 no aportan ninguna información interesante al análisis. Además, eliminar genes con bajos valores de cuentas nos evita obtener resultados falsos positivos.
##Seleccionando genes con al menos 5 cuentas por millón (cpm) en al menos 2 muestras
keep <- rowSums(cpm(counts) >= 5) >=2
table(keep)## keep
## FALSE TRUE
## 47623 13052
Es recomendable no realizar un filtrado astringente, de lo contrario se recuperan solamente genes “housekeeping”.
Tanto para edgeR como DESeq2 requerimos almacenar las cuentas y algunos metadatos en un objeto de forma de lista.
Para construir la lista de edgeR requerimos indicar cuáles son los grupos que se compararán y el número de réplicas por grupo:
## [1] "A_Control_1" "A_Control_2" "A_Control_3" "M_Control_1"
## [5] "M_Control_2" "M_Control_3" "A_Verafinib_1" "A_Verafinib_2"
## [9] "A_Verafinib_3" "M_Verafinib_1" "M_Verafinib_2" "M_Verafinib_3"
##Usar los nombres de las columnas de la matriz de cuentas para seleccionar los dos grupos experimentales
groups <- factor(sub("..$", "", colnames(counts)))
table(groups)## groups
## A_Control A_Verafinib M_Control M_Verafinib
## 3 3 3 3
Crear la lista empleando la matriz de cuentas y los grupos:
## Formal class 'DGEList' [package "edgeR"] with 1 slot
## ..@ .Data:List of 3
## .. ..$ : int [1:13052, 1:12] 338 1101 319 454 47 454 542 683 339 3096 ...
## .. .. ..- attr(*, "dimnames")=List of 2
## .. .. .. ..$ : chr [1:13052] "ENSG00000000003" "ENSG00000000419" "ENSG00000000457" "ENSG00000000460" ...
## .. .. .. ..$ : chr [1:12] "A_Control_1" "A_Control_2" "A_Control_3" "M_Control_1" ...
## .. ..$ :'data.frame': 12 obs. of 3 variables:
## .. .. ..$ group : Factor w/ 4 levels "A_Control","A_Verafinib",..: 1 1 1 3 3 3 2 2 2 4 ...
## .. .. ..$ lib.size : num [1:12] 13876857 13302895 15560526 18126099 13713631 ...
## .. .. ..$ norm.factors: num [1:12] 1 1 1 1 1 1 1 1 1 1 ...
## .. ..$ :'data.frame': 13052 obs. of 1 variable:
## .. .. ..$ genes: chr [1:13052] "ENSG00000000003" "ENSG00000000419" "ENSG00000000457" "ENSG00000000460" ...
Posteriormente los datos son normalizados empleando el método de los TMM (weighted Trimmed Mean of M-values). Para ello edgeR busca y elimina los valores atípicos (expresión absoluta de muestra como expresión relativa entre muestra) y calcula los factores de normalización. En este artículo encontrarán de manera detallada la explicación de la normalización por TMM (Evans, Hardin, and Stoebel 2018).
## group lib.size norm.factors
## A_Control_1 A_Control 13876857 0.9752559
## A_Control_2 A_Control 13302895 0.9746757
## A_Control_3 A_Control 15560526 0.9713069
## M_Control_1 M_Control 18126099 0.9686802
## M_Control_2 M_Control 13713631 0.9601787
## M_Control_3 M_Control 15051090 0.9892419
## A_Verafinib_1 A_Verafinib 13422693 1.0372465
## A_Verafinib_2 A_Verafinib 13914127 1.0222470
## A_Verafinib_3 A_Verafinib 13174976 1.0294891
## M_Verafinib_1 M_Verafinib 14548005 1.0436064
## M_Verafinib_2 M_Verafinib 14587258 1.0356998
## M_Verafinib_3 M_Verafinib 14879963 0.9976993
Es importante inspeccionar los datos y verificar que estén correctamente normalizados. Para ello, se grafica la expresión absoluta vs expresión relativa para cada muestra:
La consistencia de las réplicas la podemos verificar mediante un análisis de componentes principales (PCA) o calculando la correlación (Pearson) que existe entre las muestras:
p <- cpm(edgeRlist$counts, log = T)
p <- pca(p)
biplot(p, lab = colnames(edgeRlist$counts), pointSize = 10,
title = "PCA")cormat <- cor(cpm(edgeRlist$counts, log = T))
pheatmap(cormat, border_color = NA, main = "Correlation of replicates")Sin embargo, cuando la variación entre las réplicas es muy grande debido a “batch effects” es recomendable emplear algoritmos que permitan corregir estas variaciones sin causar grandes modificaciones a las cuentas y que dichos algoritmos consideren la distribución de los datos. Recientemente, se ha descrito una herramienta que corrige esta variación tomando en cuenta la distribución bi-nomial negativa de las cuentas (Zhang, Parmigiani, and Johnson 2020).
Para el análisis de expresión diferencial requerimos:
En esta matriz le vamos a indicar a edgeR cuáles muestras corresponden a un grupo o condición experimental. Dado que algunos experimentos pueden tener diseños muy complejos (una muestra pertenece a un tipo celular y a un tratamiento) empleamos la función interna de R, modelado de matrices:
design <- model.matrix(~0+edgeRlist$samples$group)
##El término ~0 le indica a la función no incluir una columna de intersecciones y solamente incluir tantas columnas como grupos en nuestro diseño experimental
colnames(design) <- levels(edgeRlist$samples$group)
design## A_Control A_Verafinib M_Control M_Verafinib
## 1 1 0 0 0
## 2 1 0 0 0
## 3 1 0 0 0
## 4 0 0 1 0
## 5 0 0 1 0
## 6 0 0 1 0
## 7 0 1 0 0
## 8 0 1 0 0
## 9 0 1 0 0
## 10 0 0 0 1
## 11 0 0 0 1
## 12 0 0 0 1
## attr(,"assign")
## [1] 1 1 1 1
## attr(,"contrasts")
## attr(,"contrasts")$`edgeRlist$samples$group`
## [1] "contr.treatment"
Como edgeR ajusta los datos a un modelo de distribución bi-nomial negativo (parecido a Poisson) se requiere calcular un parámetro adicional de dispersión de los datos.
edgeR calcula la dispersión a tres niveles:
edgeR en los análisis de expresión diferencial puede obviar las comparaciones entre los grupos experimentales. Sin embargo, es bueno contar con una matriz de contrastes para indicarle a edgeR cuáles son las comparaciones que queremos hacer entre los datos.
contrast <- makeContrasts(
"Cells" = "A_Verafinib - A_Control",
levels = edgeRlist$design
)
contrast## Contrasts
## Levels Cells
## A_Control -1
## A_Verafinib 1
## M_Control 0
## M_Verafinib 0
En este caso, solo tenemos dos condiciones experimentales “Células A” y “Células B”. Con esta matriz de contrastes indicamos que la comparación será buscando los genes diferencialmente expresados en B con respecto a A. Sin embargo, si tuvieramos mayor cantidad de grupos experimentales los podemos comparar, por ejemplo:
Células CT = X, Y, Z
Células Tx = P, Q, R
makeContrast("Cells" = "(P + Q + R)/3 - (X, Y, Z)/3")
El objetivo es que los coeficientes sumen 0
Los datos deben ajustarse a un modelo lineal bi-nomial negativo. Para ello se utilizará la función glmQLfit con la cual se tiene un control más robusto del “error rate”.
Se realiza la comparación para obtener los genes diferencialmente expresados entre una condición y otra. Con la función glmQLFTest, estamos probando la hipótesis nula “el valor de |lfc| del gen A es igual a 0”. Por lo tanto, los valores de p y p.adj calculados están hechos con respecto a dicha hipótesis nula.
qlf.BvsA <- glmQLFTest(fit, ##Objeto en forma de lista con los datos ajustados a un modelo bi-nomial negativo
contrast = contrast[, "Cells"]) ##Matriz de contrastes
##Obtenemos los DEG con expresión distinta de 1, valor de p menor a 0.05, corrigiendo el valor de p por el método de Benjamini-Hochberg
deg.BvsA <- decideTestsDGE(qlf.BvsA, p.value = 0.05, adjust.method = "BH", lfc = 0)
table(deg.BvsA)## deg.BvsA
## -1 0 1
## 3584 5918 3550
##Guardamos los resultados en un data.frame
DEG.BvsA <- DEGResults(qlf.BvsA)
##Añadimos una columna extra a los resultados anteriores indicando la condición de expresión respecto al lfc y al FDR
DEG.BvsA <- edgeResults(DEG.BvsA, logfc = 0, padj = 0.05)Los datos los visualizamos graficando un “Volcano plot”
Un procedimiento muy común es seleccionar aquellos genes cuya expresión difeerencial haya sido significativa y que cumplan con un valor de lfc. Por ejemplo: |lfc| > 1 (es decir genes con cuya expresión es > 2 o < 0.5 respecto al CT)
significant.1 <- DEG.BvsA %>% filter(logFC > 1 & FDR < 0.05 | logFC < -1 & FDR < 0.05)
paste(length(significant.1$genes), "is the number of significant genes")## [1] "3161 is the number of significant genes"
¡Cuidado! Filtrar los resultados de esta manera es incorrecto. Recordar que en el análisis de expresión diferencial probamos la hipótesis nula “El |lfc| del gen A es igual 0” y los valores de p y FDR corresponden a dicha prueba. El asumir que estos genes filtrados tienen |lfc| > 1 con el valor de FDR previamente calculado provoca un sesgo. Favorece genes con baja expresión y alta variabilidad, invalidando el poder estadístico de la prueba. En otras palabras, incrementamos el riesgo de captar resultados falsos positivos. Aquí encontrarán un artículo en donde se explica detalladamente las implicaciones de realizar este tipo erroneo de filtrados (Love, Huber, and Anders 2014).
Para resolver este problema, tanto edgeR como DESeq2 implementan funciones en donde se prueba la hipótesis nula “el |lfc| del gen A es distinto a x”:
qlf.BvsA.lfc1 <- glmTreat(fit,
contrast = contrast[, "Cells"],
lfc = 1)
deg.BvsA.lfc1 <- decideTestsDGE(qlf.BvsA.lfc1, p.value = 0.05, adjust.method = "BH", lfc = 1)
table(deg.BvsA.lfc1)## deg.BvsA.lfc1
## -1 0 1
## 792 11343 917
De esta manera podemos seleccionar los genes que estadísticamente tienen |lfc| > 1 y robustecer nuestros resultados.
DEG.BvsA.lfc1 <- DEGResults(qlf.BvsA.lfc1)
DEG.BvsA.lfc1 <- edgeResults(DEG.BvsA.lfc1, logfc = 1, padj = 0.05)Visualizamos estos datos en un nuevo volcano plot
significant.genes <- DEG.BvsA.lfc1 %>% filter(FDR < 0.05 & logFC > 1 | FDR < 0.05 & logFC < -1)
paste("The number of significant genes with |lfc| > 1 is", length(significant.genes$genes))## [1] "The number of significant genes with |lfc| > 1 is 1709"
Finalmente, los los genes diferencialmente expresados podemos emplearlos para crear mapas de calor “heatmaps”:
##Obtener los nombres o ids de los genes con expresión diferencial
significant.ids <- significant.genes$genes
##Crear una matriz de cuentas normalizadas por cpm empleando las cuentas que se encuentran guardadas en el objeto edgeRlist
significant.cpm <- cpm(edgeRlist$counts, log = T)
##Cortamos los genes con expresión significativa
significant.cpm <- significant.cpm[significant.ids, c(1, 2, 3, 7, 8, 9)]
##Generamos una paleta de colores
RedBlackGreen <- maPalette(low = "green", high = "red", mid = "black")
##Obtenemos el heatmap
pheatmap(significant.cpm,
border_color = NA,
color = RedBlackGreen,
show_rownames = F,
scale = "row",
angle_col = 0)Una vez obtenida la lista de genes con expresión diferencial significativa se procede a realizar análisis de representación de vías (ORA o GSEA) para darle una explicación biloógica a los resultados obtenidos.
En DESeq2 requerímos crear una tabla de metadatos (distinta a la que empleada para importar los datos) con la siguiente información:
##Crear vectores indicando el nombre de las muestras (líneas celulares) y el número de réplicas en cada condición
cell_lines <- c(rep("A", 6))
##Crear un vector indicando las dos condiciones experimentales a comparar
condition <- c(rep("control", 3), rep("tx", 3))
##En un df alojar los datos de las líneas celulares y las condiciones experimentales a las que pertenecen
meta_cells <- data.frame(cell_lines, condition)
##Especificar el nombre de las filas de acuerdo al nombre de las columnas de la matriz de cuentas
rownames(meta_cells) <- names(counts_DESeq)
meta_cells## cell_lines condition
## A_Control_1 A control
## A_Control_2 A control
## A_Control_3 A control
## A_Verafinib_1 A tx
## A_Verafinib_2 A tx
## A_Verafinib_3 A tx
Al igual que en edgeR, DESeq2 guarda o aloja los datos de las cuentas y metadatos en un objeto con clase de lista:
DDS_list <- DESeqDataSetFromMatrix(countData = counts_DESeq, ## Datos de las cuentas crudas
colData = meta_cells, ## Tabla de metadatos
design = ~condition) ## Grupos experimentales a los que pertence cada muestraEn este punto tenemos que especificarle a DESeq cuál es el grupo de referencia para realizar las comparaciones:
Para filtrar los genes con expresión baja, DESeq emplea una función interna de normalización para que los datos de las cuentas sean comparables entre ellos:
##Usamos los mismos criterios que en edgeR
filtered <- rowSums(counts(DDS_list) >= 5) >= 2
table(filtered)## filtered
## FALSE TRUE
## 41753 18922
Las muestras son normalizadas empleando la función DESeq. Esta función realiza tres pasos en un solo comando:
## estimating size factors
## estimating dispersions
## gene-wise dispersion estimates
## mean-dispersion relationship
## final dispersion estimates
## fitting model and testing
Finalmente, obtenemos los genes diferencialmente expresados. En este caso vamos a probar la hipótesis nula “El lfc del gen A es igual a 0”:
##
## out of 18922 with nonzero total read count
## adjusted p-value < 0.05
## LFC > 0 (up) : 4029, 21%
## LFC < 0 (down) : 4162, 22%
## outliers [1] : 8, 0.042%
## low counts [2] : 0, 0%
## (mean count < 2)
## [1] see 'cooksCutoff' argument of ?results
## [2] see 'independentFiltering' argument of ?results
Visualizamos los datos en un volcano plot:
##convertimos los resultados a un df
DEG.VerafinibvsCT <- as.data.frame(DEG.VerafinibvsCT) %>%
rownames_to_column(var = "ensgene_id")
DEG.VerafinibvsCT <- DESEqResults(DEG.VerafinibvsCT, logfc = 0, FDR = 0.05)
volcano_DESeq2(DEG.VerafinibvsCT, lfc = 0, FDR = 0.05)Y de manera similar a edgeR, generamos un heatmap con los genes diferencialmente expresados:
significant.genes.2 <- DEG.VerafinibvsCT %>% filter(log2FoldChange > 0 & padj < 0.05 | log2FoldChange < 0 & padj < 0.05)
significant.ids.2 <- significant.genes.2$ensgene_id
significant.counts <- counts(DDS_list, normalized = T)
significant.counts <- significant.counts[significant.ids.2, ]
pheatmap(significant.counts,
border_color = NA,
color = RedBlackGreen,
show_rownames = F,
scale = "row",
angle_col = 0)La lista de genes diferencialmente expresados obtenida por edgeR o DESeq2 no aportan información sobre los procesos biológicos que pueden verse afectados en las células tratadas comparado con el control. Para ello, se requiere realizar análisis que le den un contexto biológico a las listas obtenidas. Por ejemplo los análisis de vías son muy útiles para contextualizar los resultados:
En este artículo encontrarán la información detallada de los análisis de vías, sus fortalezas, debilidades y alcances de cada método (García-Campos, Espinal-Enríquez, and Hernández-Lemus 2015).
Dobin, Alexander, Carrie A. Davis, Felix Schlesinger, Jorg Drenkow, Chris Zaleski, Sonali Jha, Philippe Batut, Mark Chaisson, and Thomas R. Gingeras. 2013. “STAR: Ultrafast Universal RNA-Seq Aligner.” Bioinformatics 29 (1): 15–21. https://doi.org/10.1093/bioinformatics/bts635.
Evans, Ciaran, Johanna Hardin, and Daniel M. Stoebel. 2018. “Selecting Between-Sample RNA-Seq Normalization Methods from the Perspective of Their Assumptions.” Briefings in Bioinformatics 19 (5): 776–92. https://doi.org/10.1093/bib/bbx008.
García-Campos, Miguel A., Jesús Espinal-Enríquez, and Enrique Hernández-Lemus. 2015. “Pathway Analysis: State of the Art.” Frontiers in Physiology 6: 383. https://doi.org/10.3389/fphys.2015.00383.
Li, Bo, and Colin N. Dewey. 2011. “RSEM: Accurate Transcript Quantification from RNA-Seq Data with or Without a Reference Genome.” BMC Bioinformatics 12 (August): 323. https://doi.org/10.1186/1471-2105-12-323.
Love, Michael I, Wolfgang Huber, and Simon Anders. 2014. “Moderated Estimation of Fold Change and Dispersion for RNA-Seq Data with DESeq2.” Genome Biology 15 (12): 550. https://doi.org/10.1186/s13059-014-0550-8.
Zhang, Yuqing, Giovanni Parmigiani, and W Evan Johnson. 2020. “ComBat-Seq: Batch Effect Adjustment for RNA-Seq Count Data.” NAR Genomics and Bioinformatics 2 (3): lqaa078. https://doi.org/10.1093/nargab/lqaa078.